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Various instabilities have been proposed as candidates to prompt the condensation of 
■^j- ■ giant, star-forming cloud complexes from the diffuse interstellar medium. Here, we use 

three-dimensional ideal MHD simulations to investigate nonlinear development of the 
qq . Parker, magneto- Jeans (MJI), and swing mechanisms in galactic disk models. The disk 

models are local, isothermal, and begin from a vertically-stratified magnetohydrostatic 
equilibrium state with both gaseous and stellar gravity. We allow for a range of surface 
densities and rotational shear profiles, as well as unmagnetized control models. We first 
Qh' construct axisymmetric equilibria and examine their stability. Finite disk thickness 

reduces the critical Toomre stability parameter below unity; we find Q c ~ 0.75, 0.72, 
and 0.57 for zero, sub-equipartition, and equipartition magnetic field cases, respectively. 
We then pursue fully three-dimensional models. In non-self-gravitating cases, the peak 
. mid-disk density enhancement from the "pure" Parker instability is less a factor of two. 

The dominant growing modes have radial wavelengths \ x comparable to the disk scale 
height H, much shorter than the azimuthal wavelength (A^ ~ 10 — 20H). Shearing disks, 
being more favorable to midplane-symmetric modes, have somewhat different late-time 
magnetic field profiles from nonshearing disks, but otherwise saturated states are similar. 
Late-time velocity fluctuations at 10% of the sound speed persist, but no characteristic 
structural signatures of Parker modes remain in the new quasi-static equilibria. In 
self-gravitating cases, the development of density structure is qualitatively similar to 
our previous results from thin-disk simulations. The Parker instability, although it 
may help seed structure or tip the balance under marginal conditions, appears to play a 
secondary role - not affecting, for example, the sizes or spacings of the bound structures 
that form. In shearing disks with Q less than a threshold level ~ 1, swing amplification 
can produce bound clouds of a few times the local Jeans mass. The most powerful 
cloud-condensing mechanism, requiring low-shear conditions as occur in spiral arms or 
galactic centers, appears to be the MJI. In thick disks, the MJI occurs for \ y ^ 2ttH. 
Our simulations show that condensations of a local Jeans mass ( < 3 x 10 7 Mq) grow 
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very rapidly, supporting the idea that M JI is at least partly responsible for the formation 
of bound cloud complexes in spiral galaxies. 

Subject headings: galaxies: ISM — galaxies: kinematics and dynamics — galaxies: 
structure — instabilities — ISM: kinematics and dynamics — ISM: magnetic fields - 
MHD — stars: formation 

1. Introduction 

Nearly half of the interstellar medium (ISM) in the inner part of the Milky Way is estimated 
to be in the molecular component (e.g., Dame 1993), with the largest portion of the molecular 
ISM in giant molecular clouds or cloud complexes (GMCs) of masses ~ 10 4 — 10 6 M (Dame et al. 
1987; Solomon et al. 1987), where most of star formation occurs (e.g., Williams & McKee 1997). 
GMCs are generally turbulent, self-gravitating, and magnetized (e.g., Blitz 1993; Crutcher 1999). 
In external galaxies, they often appear in clusters, forming giant molecular associations (GMAs) 
(e.g., Grabelsky et al. 1987; Vogel et al. 1988; Rand & Kulkarni 1990; Rand 1993; Sakamoto et al. 
1999), or within the boundaries of H I superclouds (e.g., Elmegreen 1995). These giant clouds in 
spiral galaxies are mostly associated with spiral arms, but they sometimes appear in the interarm 
regions (e.g., Solomon et al. 1985; Rand 1993; Kenney 1997; Heyer &: Terebey 1998). By regulating 
the rate and mode of star formation in disk galaxies, GMCs play a fundamental role in controlling 
galactic evolution and may be a key to understanding the nature of the Hubble sequence. But how 
GMCs form is still not clearly understood. 

Traditionally, proposed mechanisms for giant cloud formation fall into two categories: stochas- 
tic coagulation of smaller clouds (following Oort 1954; e.g., Kwan 1979) or collective effects involving 
instabilities (see e.g., Elmegreen 1996). Existing work suggests that collisional agglomeration may 
proceed too slowly (e.g., Scoville & Hersh 1979; Elmegreen 1990, 1995). Collisions in most cases 
result in disruption rather than merger (e.g., Lattanzio et al. 1985; Kim et al. 1999; Klein, Woods, 
& McKee 2001); the available mass in smaller clouds is not enough to form GMCs through colli- 
sional build-up (Heyer & Terebey 1998; Blitz & Williams 1999); collisional agglomeration would 
lead to magnetically supercritical clouds and active star formation before GMC-scale masses are 
reached (Ostriker, Gammie, & Stone 1999; Ostriker, Stone, & Gammie 2001); GMC chemistry may 
be inconsistent with an extended build-up phase (e.g., van Dishoek et al. 1993). All of these and 
other difficulties point to the need for collective effects (Blitz &: Shu 1980). This second means 
of building GMCs and GMAs involves large-scale dynamical instabilities that coherently amass 
material over scales large compared to the sizes of diffuse atomic clouds (see e.g., Elmegreen 1995). 
The small-scale structure of the medium plays a secondary role, with the modal gathering of ma- 
terial operating collectively on the whole distribution. Coherently-growing unstable modes that 
produce condensations depend on various effects, including self-gravity, the Coriolis force, sheared 
azimuthal motion, magnetic tension, and magnetic buoyancy, all of which have been shown to 
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interact in unstable fashion. 

The two most well-known condensation instabilities are the Parker instability (Parker 1966), 
in which buoyancy causes dominantly in-plane magnetic fields to buckle, with matter collecting 
in valleys; and the swing amplifier, in which epicyclic motion conspires with galactic shear to en- 
able self-gravitating enhancement of overdense regions (Goldreich & Lynden-Bell 1965b; Toomre 
1981). Another condensation instability that has been shown to be efficient in low-shear, magne- 
tized regions is the magneto- Jeans instability (MJI), in which magnetic tension of in-plane fields 
counters the stabilizing effect of the Coriolis force, allowing self-gravitating contraction of over- 
dense regions (Elmegreen 1987; Kim & Ostriker 2001). Because the growth rates of Parker and 
MJI modes increase with magnetic field strength and density, the compression of the ISM induced 
by passage through stellar spiral arms may be important in triggering GMC and GMA formation 
(e.g., Mouschovias, Shu, & Woodward 1974; Blitz & Shu 1980). Observed distributions of H II 
regions and OB star complexes in a "beads on a string" pattern along spiral arms may be the 
direct consequence of cooperation among the Parker instability, Jeans instability, and spiral arm 
potentials (e.g., Elmegreen & Elmegreen 1983). 

In Kim & Ostriker (2001, hereafter Paper I), we presented studies of detailed nonlinear evolu- 
tion of self-gravitating modes in shearing, razor-thin galactic disks. We showed that formation of 
bound condensations via the swing instability is subject to "threshold" behavior, with models hav- 
ing sufficiently small values of the Toomre Q parameter (see eq. [8] below) producing bound clouds, 
and models with larger values of Q remaining nonlinearly stable. We found that the threshold value 
of Q is in the range 1.2 — 1.4 (for magnetic field strengths from zero up to thermal-equipartition 
values), consistent with an apparent mean threshold Q t h = 1-4 for star formation evident from ob- 
servational studies (Kennicutt 1989; Martin & Kennicutt 2001). In Paper I, we also elucidated the 
role of MJI in forming bound condensations by showing that in low-shear regions where the swing 
amplifier cannot operate, the mechanism that leads to growth of perturbations is quite different 
from that in high-shear regions, and that the presence of magnetic fields is crucial to the instability. 
In the absence of magnetic fields, specific vorticity is conserved; contraction of vorticity-conserving 
masses leads to increased Coriolis forces that exceed the increase in self-gravity such that further 
contraction is prevented (in the absence of a shearing background). In magnetized systems, how- 
ever, tension forces can transfer vorticity from one region to another, such that condensation driven 
by self-gravity can proceed. 

Effects of spiral arm gravity on the growth of giant clouds were explicitly addressed by Kim & 
Ostriker (2002, hereafter Paper II) using two-dimensional models. Paper II showed that for a spiral 



1 Thermal instabilities involving changes in the microphysical kinetic temperature are important in forming cold 
clouds on small scales (Field 1965; see also e.g., Hennebelle & Perault 1999 and Burkert & Lin 2000 for recent 
simulations) ; if dissipation of kinetic energy in collisions of cloudlets acts effectively like microscopic cooling for the 
turbulent "cloud fluid", then analogous mechanisms may also be important on large scales (e.g., Elmegreen 1989, 
1991). 
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arm inside corotation, the enhanced surface density and shear gradients in the arm region promote 
the growth of gaseous spurs. These spurs jut nearly perpendicularly from the outer (downstream) 
side of the arm, and trail in the same sense as the arm itself; this is consistent with observed 
spur-like dust lanes in M51 (Scoville et al. 2001). We argued that the formation of spurs can 
be understood as a natural consequence of the MJI within spiral arms, where compression leads 
to a reduced (or negative) shear gradient that allows these self-gravitating modes to grow over 
an extended time. Paper II showed that bound clouds of typical mass ~ 4 x 10 6 M Q can form 
from spur fragmentation; such condensations potentially represent structures that would initially 
be observed as GMAs, and subsequently evolve to become bright OB associations. 

Because of their restriction to razor-thin geometry, our previous two-dimensional models could 
not capture the potential consequences of the Parker instability on cloud formation. While previous 
work on the Parker instability has studied the weak-shear limit (see linear-theory analysis of Shu 
1974, Foglizzo Sz Tagger 1994, 1995) and the case of solid-body rotation (see linear theory of 
Zweibel & Kulsrud 1975 and simulations by Chou et al. 2000 and Kim, Ryu, & Jones 2001), and 
indeed shown that rotation and shear tend to be stabilizing, no previous study focused primarily 
on Parker modes has incorporated shear at a realistic "average" galactic level. Although the Parker 
instability has been proposed to prompt the formation of molecular clouds (e.g., Mouschovias, Shu, 
&; Woodward 1974; Blitz & Shu 1980), existing numerical simulations indicate that the Parker 
instability alone is unable to produce structures of high enough surface density to represent GMCs 
(Basu, Mouschovias, & Paleologou 1997; J. Kim et al. 1998, 2000; Santillan et al. 2000; Kim, Ryu, 
& Jones 2001). 2 

To study the coupling of the Parker instability with self-gravitating modes, Elmegreen (1982a,b) 
performed linear stability analyses for combined Parker-Jeans modes and concluded that clouds of 
mass ~ 10 5 — 10 6 M Q form in the density regime for which self-gravity significantly modifies the 
Parker instability; this regime of densities is expected to be found in spiral arm regions only. 
Hanawa, Nakamura, & Nakano (1992) included solid-body rotation in a linear-theory analysis of 
the Parker-Jeans instability, and concluded that the transfer of angular momentum out of growing 
condensations may be a limiting factor. This transfer of angular momentum may be accomplished 
by magnetic stresses much as it occurs in the two-dimensional MJI. A nonlinear simulation of the 
Parker-Jeans instability in a disk with solid-body rotation was among the models performed by 
Chou et al. (2000), who found that filamentary structures tend to form perpendicular to the mean 
magnetic field. 

To our knowledge, there have previously been neither linear-theory nor numerical studies that 
incorporate all of the aforementioned elements that can contribute to three-dimensional instabilities 
in galactic disks: buoyant (horizontal) magnetic fields, self-gravity, rotation, and shear. In this 
paper, we investigate nonlinear evolution of galactic disks allowing for all these fundamental physical 



2 In addition, a random component of galactic magnetic fields whose strength is almost comparable to the mean 
uniform component is shown to suppress the Parker instability significantly (Kim & Ryu 2001). 
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agents. This work extends our previous studies (Paper I) of structure formation via gravitational 
instability by including the effects of the Parker instability. Our primary objectives are to study 
how the Parker instability develops in the presence of strong shear, and to compare the eventual 
evolution of disks under the various instabilities separately or in combination. We can thus assess 
the potential for gravitationally-bound cloud formation from a variety of dynamical mechanisms. 
Since vertically-integrated, razor-thin disks are known to overestimate self-gravity at the midplane 
(e.g., Toomre 1964), we will also study analytically and numerically how the stability of three- 
dimensional disks to in-plane self-gravitating modes changes as their thickness varies. 

In our model, disks are local (see §2.1), isothermal in both space and time, and vertically 
stratified. Initial magnetic fields are plane-parallel, pointing in the azimuthal direction. To model 
realistic conditions in disk galaxies, we include both external gravity arising from stars (as a time- 
independent potential) and gaseous self-gravity, which together counterbalance thermal and mag- 
netic pressure gradients in an unperturbed state. We do not, however, consider in this paper other 
features such as stellar spiral arms and effects of cosmic rays that may enhance the magnetic buoy- 
ancy force. We begin by solving the governing equations to obtain initial static equilibrium density 
distributions, and then explore their gravitational stability to axisymmetric perturbations. We 
next select several parameter sets that represent "average disk" or "spiral arm" regions. After ap- 
plying small-amplitude perturbations to the corresponding density profiles, we follow the nonlinear 
dynamical development with three-dimensional direct numerical simulations. 

The paper is organized as follows: In §2, we describe the basic MHD equations in the local 
"shearing-sheet" approximation, the computational methods we use, and our model parameters. In 
§3, we present initial magnetohydrostatic equilibria in the presence of both self-gravity and external 
gravity, and analytically examine their axisymmetric stability by making simplifying assumptions. 
In §4, we present the results of three-dimensional simulations for the Parker instability without 
self-gravity (§4.1), generalized (Parker /swing) gravitational instability in shearing regions (§4.2), 
and MJI in nonshearing regions (§4.3). We summarize our results and discuss the implications of 
present work for galactic cloud formation in §5. 

2. Numerical Methods and Model Parameters 

2.1. Basic Equations and Numerical Methods 

In this paper, we study nonlinear evolution of vertically stratified, differentially rotating, self- 
gravitating galactic gas disks with threaded magnetic fields in the azimuthal direction. We solve 
fully three-dimensional, compressible, ideal MHD equations in a local Cartesian reference frame 
whose center lies at a galactocentric radius Ro and orbits the galaxy with a fixed angular velocity 
SI = Q(Rq). In this local frame, radial, azimuthal, and vertical coordinates are represented by 
x = R — Ro, y = Ro(<fi — Qt), and z, respectively, and terms associated with the curvature effect are 
neglected (Goldreich & Lynden-Bell 1965b; Julian & Toomre 1966). The local- frame equilibrium 
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background velocity relative to the center of the box atx = y = z = 0is then reduced to 
v o = — Q^xy, where 

din ft 
q= — 



dlnP 



(1) 

Ro 



is the local dimensionless shear rate. In terms of q, the local epicyclic frequency k is given by 
k 2 = (4 — 2q)Q 2 . The resulting "shearing-sheet" equations expanded in the local frame are 

^ + V-(pv)=0, (2) 

-J + v • Vv = — VP + - — (V x B) x B + 2qn 2 xx - 2ft x v - V$ s + g ext , (3) 
at p 4-Kp 

<9B 

— = V x (v x B), (4) 
V 2 $ s = AirGp, (5) 

and 

P = c^, (6) 

(cf, Hawley, Gammie, & Balbus 1995; Stone et al. 1996). Here, v is the velocity in the local frame, 
p is the mass density, P is the thermal pressure, B is the magnetic field, <I> S is the self-gravitational 
potential, g cx t is the external gravitational acceleration, G is the gravitational constant, and c s 
is the isothermal sound speed. In equation (3), the third and fourth terms on the right hand 
side represent the tidal and Coriolis forces, respectively, which balance each other in the initial 
equilibrium with v = vo- As expressed by equation (6), for simplicity we adopt an isothermal 
equation of state for all the simulations presented in this paper. 

An initial vertical equilibrium is defined by the balance between the total (thermal plus mag- 
netic) pressure gradient and the total (self plus external) gravitational force. For g ext representing 
the external gravity from stars, we adopt the simple form 

gext = -T7-Z, (7) 
-T2* 

where H* is the scale height of stars and is a free parameter. Since most of gas is distributed 
within \z\ < P* and the small- 1 2 1 part of the realistic gravity deduced from stellar tracer populations 
can be described by a linear function of height (e.g., Bahcall 1986; Kuijken & Gilmore 1989), g cx t 
given in equation (7) is a fairly good approximation in studying dynamics of the gaseous component. 
Although a linear gravity law is known to strengthen the Parker instability slightly compared to 
the case with a more realistic gravity law, the growth rate of the most unstable mode differs by 
only a factor of 1.2 (Kim & Hong 1998). As we shall show in §3, the total gravity that drives the 
Parker instability in our models is similar to the external gravity adopted by other authors (Giz & 
Shu 1993; Kim & Hong 1998). 
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The numerical solutions to the dynamical equations (2) — (7) are obtained using a modified 
version of the ZEUS code originally developed by Stone & Norman ( 1992a, b). ZEUS is a time- 
explicit, operator-split, finite-difference method for solving the MHD equations on a staggered mesh, 
and has been demonstrated to provide robust solutions to a wide variety of astrophysical problems. 
ZEUS employs "constrained transport" to guarantee that V • B = within machine precision, 
and the "method of characteristics" for accurate propagation of Alfvenic disturbances (Evans & 
Hawley 1988; Hawley & Stone 1995). For less diffusive advection of hydrodynamic variables, we 
apply a velocity decomposition method which treats the contribution from the background shearing 
parts as source terms, while solving the transport step using only the perturbed parts. In order to 
minimize the errors arising from discontinuities in the flow characteristics across the x-boundaries, 
the ghost zones adjoining the x-boundaries are kept active (Paper I). In the x and y directions, 
we implement shearing box boundary conditions in which the x-boundaries are shearing-periodic 
and the y-boundaries are perfectly periodic (Hawley, Gammie, & Balbus 1995). In the z direction, 
we adopt free (outflow) boundary conditions which set all the variables in the boundary zones 
equal to the corresponding values in the adjacent active zones (Stone & Norman 1992a). These 
boundary conditions permit outgoing waves and thus physical quantities such as mass, momentum, 
and magnetic flux can flow out of the simulation domain. Since density is very low near the vertical 
boundaries, however, simulation results generally show that the total mass is conserved within 1%. 

For computation of the self-gravitational potential, we employ a hybrid technique that combines 
the Green function method to compute the gravitational kernels (Miyama, Narita, & Hayashi 
1987) and the Fourier transform method in shearing-sheet coordinates (Gammie 2001). Our high- 
resolution models have a 256 x 128 x 128 grid for (x, y, z); we find that the Poisson solver typically 
takes about 48% of the total cpu time. We have checked our implementation of the ZEUS code 
on various test problems including two-dimensional gravitational instability and propagation of 
gravito-MHD waves in three dimensions. 



2.2. Model Parameters 

Our model disks are highly-idealized, smoothed versions of what is physically a much more 
complex, multiphase interstellar medium (e.g., Field, Goldsmith, & Habing 1969; Cox & Smith 1974; 
McKee &: Ostriker 1977; Heiles 2001). We adopt an isothermal condition in both space and time, 
with an effective isothermal speed of sound c s = 7 km s -1 , corresponding to a mean Galactic thermal 
pressure P/k ~ 2000 — 4000 Kcm~ 3 (Heiles 2001) and mean midplane density «h ~ 0.6 cm -3 
(Dickey & Lockman 1990); external galaxies have similar gaseous velocity dispersions. One of 
the parameters in our simulations is the rotation rate 0. When reporting dimensional results, we 
normalize to the solar neighborhood value O = 26kms~ 1 kpc~ 1 (Binney & Tremaine 1987). For a 
flat rotation curve, q ~ 1 and thus k ~ 2 1 / 2 0; we allow for varying q to represent both "average" 
regions and reduced-shear (spiral arm) regions. 

Our initial disk models are vertically stratified, with initial density distribution po(z). The 
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corresponding surface density is given by So = po{z)dz. To characterize So in a dimensionless 
fashion, we use the Toomre Q stability parameter 



For an infinitesimally-thin disk, the critical value for axisymmetric gravitational stability is Q c = 1 
(Toomre 1964; see also e.g., Binney &: Tremaine 1987; Shu 1992). We shall show in §3 that the finite 
thickness of disks reduces Q c (cf. Goldreich & Lynden-Bell 1965a). For local gas surface density of 
So = 13 M pc~ 2 comparable to the solar neighborhood (cf, Paper I; Kuijken & Gilmore 1989), 
Q « 1.5; we adopt this value of Q for our "average disk" fiducial value. 

The initial magnetic field in our models points in the azimuthal (y) direction and varies with 
height as Bo = Bo(z)y; its magnitude is parameterized by 

0=4-> ( 9 ) 
V A 

where the Alfven speed is va = Bo/y/Anpo. As an initial condition, we fix (3 = 1 spatially for all 
the MHD models in this paper, so that initially the Alfven speed is the same as the isothermal 
sound speed everywhere. This implies Bq(z) oc \J po{z). 

We explore the nonlinear evolution of disks under both average disk and other galactic con- 
ditions by changing So and q. In order to model a portion of a gas disk inside a spiral arm, we 
use So = 65 M Q pc~ 2 , which is five times larger than the average disk value. Since k oc S 1 / 2 and 
Q oc S^ 1 / 2 for isothermal gas from the constraint of potential vorticity conservation (Balbus & 
Cowie 1985; Paper II), this higher surface density corresponds to Q ~ 0.7, which we adopt as a 
fiducial "spiral arm" value. We allow the shear parameter q to vary from model to model, but fix 
it spatially in a given model. Note that this arbitrary implementation of background shear is not 
consistent with realistic shear that varies rapidly across the spiral arm (Paper II). Nevertheless, 
our high surface density models allow us to focus on the effects of enhanced self-gravity inside the 
spiral arm without the additional complications of in-plane variations of the background. 

We parameterize the magnitude of the external gravity by 



which measures the relative strength of self- (gaseous) gravity and external (stellar) gravity in the 
equilibrium state (cf, Elmegreen 1982a). Here, p^ denotes the midplane density in the absence of 
both magnetic field and external gravity (i.e., (3 = oo and sq = oo; see §3). If the stellar distribution 
approximates a vertical self-gravitating equilibrium with vertical velocity dispersion <7* iZ , one can 
rewrite equation (10) as 



iirGpooH* 



(10) 



s = 



9* 




(11) 
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where denotes the stellar surface density. 3 Notice that so is now expressed in terms of observable 
quantities. Since ~ 35 M pc -2 and <7* )2 ~ 20 kms -1 (e.g., Kuijken & Gilmore 1989; Holmberg 
&: Flynn 2000), so ~ 1 is suitable for "average disk" conditions, while sq ~ 25 would describe a 
spiral arm region (taking the stellar spiral perturbation as weak). 

Our simulation domain is a rectangular parallelepiped of size L x x L y x L z , with a vertical 
range of \z\ < L z /2 = 4H, where H is the gaseous scale height (see eq. [14] below). The horizontal 
dimensions of the box are L x = L y = 25H or 17H for so = 1 or so = 25 models, respectively. 
This size corresponds to the most unstable wavelength of the axisymmetric gravitational instability 
in each model with (3=1 (see §3.2). Conversion of H into a physical value requires solving 
equations governing magnetohydrostatic equilibria, which will be presented in the next section. 
Using c s /£l = 270pc (c s /7kms _1 )(r2/26kms _1 kpc -1 ) -1 and borrowing the results from Figure 
1, we obtain for high surface-density (spiral arm) regions (so = 25), H ~ 50 pc for (3 = oo and 
H rs 80 pc for (3 = 1 models. Under average disk conditions (so = 1), H « 160 pc for [3 = oo and 
H « 210 pc for (3 = 1 cases. The former value for an unmagnetized model is consistent with the 
usual inner-galaxy scale height of ~130-170 pc derived from the hydrostatic force balance between 
external gravity and turbulent pressure in the average disk condition (see e.g., Boulares & Cox 
1990; Lockman & Gehman 1991; J. Kim et al. 2000). On the other hand, the latter value for a 
(3=1 model is similar to the direct estimate ~220 pc for the H I scale height obtained by tangent 
point observations near the solar circle (Malhotra 1995), beyond which the H I disk is known to 
flare dramatically (Kulkarni, Blitz, & Heiles 1982). 

The larger initial scale height of magnetized models reflect the potential importance of mag- 
netic pressure in supporting a galactic disk against gravity. If magnetic support is significant but 
is neglected in observationally-based determinations of the scale height from surface density and 
velocity dispersion, then H could be underestimated. We note, however, that the dynamical rear- 
rangement of magnetic flux in our simulations always results in a reduction of the scale height over 
time, so that the initial values of H in our (3 = 1 models is an upper limit for a real disk of high 
magnetization. 

The total mass M to t = ^oL x L y contained in the simulation box corresponds to ~ 3.6 x 
10 8 M Q (#/210pc) 3 for s = 1, (3 = 1 models, ~ 1.2 x 10 8 M Q (#/80pc) 3 for s = 25, (3=1 models, 
and ~ 0.5 x 10 8 M Q (#/50pc) 3 for s = 25, (3 = oo models. 

For all the simulations reported in this paper, we apply initial perturbations only to the density. 
We first generate a Gaussian random variable 5 using the usual Box-Muller algorithm (cf, Press et 
al. 1992) and then construct its power spectrum in the Fourier space such that (|5a;| 2 ) oc fc~ n / 3 for 
1 < kL x /2ir < 32 and zero for 32 < kL x /2ir. We normalize by fixing the standard deviation of S in 
real space to be 1%. The final form of density perturbations is given by the product of 8 and the 
background density profile po. 



3 For an application to a central-mass dominated Keplerian disk, g exi = — Q 2 z. Then, so = 47rGp 00 /n 2 = 2Q 
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Finally, we remark on a few dynamical timescales of note. The rotation time is denned by 
trot = 1/^ = 3.8 x 10 7 yrs (r2/26kms~ 1 kpc -1 ) -1 . The shearing time is t s h = which is 

equal to t ro t for a flat rotation curve (q = 1). Related to t s h is the swing amplification time 
scale, which is approximately t sw ~ (4 — 5)t s h (cf, Toomre 1981; Paper I). The orbital period 
is t or b = 2ir/Q = 2.4 x 10 8 yrs (17/26 km s -1 kpc -1 ) -1 , which is the time unit we adopt for our 
presentation. The typical growth time of the Parker instability is given by the Alfven crossing time 
over the scale height 

tp „ kor ^= M x 10' yrs (^) (^) , (12) 
while a characteristic gravitational contraction time is 

(Paper I). Note that the Parker instability involves vertical motions of gas, while fluid motions in 
the gravitational instability are mostly parallel to the galactic plane at least in the linear regime. 



3. Initial Equilibria and Axisymmetric Stability 

3.1. Initial Magnetohydrostatic Equilibria 

In the absence of perturbations, our isothermal disk is in vertical magnetohydrostatic equi- 
librium with a density profile po(z)- Thermal and magnetic pressure forces balance against the 
combination of gaseous self-gravity and external gravity. Stronger gravity leads to a thinner disk. 
In this section, we study the change of disk thickness as sq varies. 

For simplicity, let us define poo be the midplane value of po and 

p (z)dz, (14) 



* = -/ 

Poo Jo 



be the gaseous scale height for general cases, while p^ and are the counterparts for purely 
self-gravitating, unmagnetized disks (i.e., so = oo and (3 = oo). From the definitions, it follows for 
fixed So that H = So/(2poo) = -HooAx>/poo an d using the Spitzer (1942) isothermal disk solution 



(2vrGp 00 ) 1 /2 vrGS VTkms- 1 / V 13M oP c 

Arranging equations (3) and (5) - (7) for a static equilibrium in the vertical direction, we are left 
with 

If 1\ d 2 lnp 1 

2\} + w)-d?- = -^r p -^ (16) 
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where the dimensionless height ( = z/H^ and the dimensionless density p{z) = po(z)/poo- Equa- 
tion (16) is to be integrated subject to the constraint (14). 

In the two limiting cases of so, equations (14) and (16) yield analytic solutions. For s S> 1 
corresponding to a strongly self-gravitating disk and negligible external gravity, we obtain p = 
sech 2 (z/H) along with H/H^ = 1 + 1/2/3 (e.g., Spitzer 1942; Elmegreen & Elmegreen 1978), while 
for so <C 1 corresponding to weak self-gravity compared to external gravity, p = exp(— ttz 2 /AH 2 ) 
along with HjH^ = (7rso) 1 ^ 2 (l + l/2/3) 1 / 2 /2. Note that the latter Gaussian function is the usual 
solution for density in an equilibrium under linear gravity, with AH 2 /ir = 2c 2 H if /g Je . 

We solve equations (14) and (16) iteratively: for general so, we arbitrarily choose an initial 
H/Hoo and numerically integrate equation (16) to find /5(C)- If the calculated p{C) does not satisfy 
the constraint (14), we modify H/H^ and iterate until convergence is attained. The resulting 
H/Hco as a function of sq is shown in Figure 1 where the solid line corresponds to the unmagnetized 
case (f3 = oo), while dotted and dashed lines are for sub-equipartition (f3 = 10) and equipartition 
{(3 = 1) magnetic field cases, respectively. To compare H with a length scale independent of the 
gas disk's mass, Figure 1 also plots HQ/c s = (H/H oc )(QQ,/k) as a function of sq. One may think 
of increasing sq as either decreasing the stellar content (£*) while holding the gas content (So) 
fixed, or increasing So for a given value of S*. Since is dependent on So, the first meaning 
of so is helpful in interpreting the H/H^ curves: moving toward larger (smaller) so, a gas disk 
thickens (thins) as the stellar content relative to gas decreases (increases). The magnetic pressure 
support has a larger impact at higher so- On the other hand, since £l/c s is independent of so and 
assumed to be constant at a given radius in a galaxy, the second meaning of so helps interpret 
the H£l/c s curves: as gas moves from interarm regions into spiral arms where density and so are 
higher, the gas layer becomes thinner in response to its own enhanced self-gravity. For instance, 
H(sq = 1)/H(sq = 25) w 3 when j3 = oo. Note that to take into account the effects of compressed 
magnetic fields, (3 should change correspondingly in spiral arms. 

As mentioned before, the effective driving force, g e s, of the Parker instability in our models is 
the sum of the external gravity and the self-gravity of a disk, 

9cS = ffcxt - -^7- (17) 

In Figure 2 we plot g c g as solid lines together with contributions of gaseous self-gravity (dotted 
lines) and external stellar gravity (short-dashed lines). For average disk conditions (Q ~ 1.5 and 
So ~ 1), for \z\/H ^, 1, external gravity is as important as self-gravity in establishing an equilibrium 
(Fig. 2a). For comparison, we also plot the external gravity models adopted by Giz & Shu (1993) 
and Kim k, Hong (1998), which are in good agreement with our total average disk gravity model. 
On the other hand, strongly-compressed spiral arm regions are entirely dominated by self-gravity 
(Fig. 26). 
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3.2. Axisymmetric Stability 



Dynamics of galactic gaseous disks have often been studied by using height-integrated equa- 
tions, accurate when the vertical height of the disk is much smaller than the physical scale of interest 
(e.g, Binney & Tremaine 1987; Shu 1992). One of the drawbacks of this thin-disk approximation is 
that it overestimates self-gravity for horizontal scales approaching the disk thickness (e.g., Toomre 
1964). As mentioned earlier, the thin-disk approximation leads to Q c = 1 for axisymmetric gravita- 
tional instability. Goldreich & Lynden-Bell (1965a) studied the stability of uniformly rotating disks 
of finite thickness, and showed that Q c = 0.676 for a purely self-gravitating (so = oo), isothermal 
disk (see also Gammie 2001). In this section, we explore the axisymmetric stability of isothermal 
disks under both self-gravity and external gravity. 

For analytic simplicity, we make an approximation that fluid motions are all lateral, confined 
to the x-y plane. This ignores the small vertical oscillations of the disk as well as the potential 
effects of the Parker instability (which would be nonaxisymmetric anyway). We perform a standard 
linear stability analysis for axisymmetric modes, assuming perturbations are a linear superposition 
of Fourier modes, exp(iu;i — ikx), of frequency to and wavenumber k. For the gravitational force 
due to a perturbed density distribution pi = pi(z)e~ ikx with k ^ 0, the Poisson equation (5) yields 
an approximate solution d&i/dx = Am sgn(k)Gpi(0)He~ tkx /(l + \k\H) at the disk midplane. This 
expression for the perturbed force is exact for an exponential density distribution pi(z) oc e~ z l H 
(e.g., Elmegreen 1987), and asymptotically approaches the exact solutions for arbitrary density 
profiles when \k\H — ► or oo. For the vertical density profiles that we found in the previous 
section, the largest fractional difference between the true and approximate values for the perturbed 
radial gravity is only ~ 15%, which occurs at \k\H ~ 1. The procedure for driving the dispersion 
relation is quite straightforward and we write down only the result, 



For short-wavelength perturbations kH » 1, equation (18) recovers the dispersion relation of waves 
in a rotating medium with uniform density poo in three dimensions (Chandrasekhar 1961). For long- 
wavelength perturbations kH <C 1, equation (18) reduces to the dispersion relation in a razor-thin 
disk with surface density So = 2pooH. 

Using the H — so relationship that we found in the previous subsection, we compute fc max 
which minimizes the right hand side of equation (18) for given so an d So (and thus maximizes the 
growth rate). The results are shown in Figure 3a. The most unstable wavelength for axisymmetric 
perturbations is typically A max = 27r/fc max ~ 20H for average disk conditions (i.e., so = 1) and 
decreases slightly as gaseous self-gravity becomes more important (larger sq). For fixed gas content, 
the /cmaxifoo curves show that relatively weaker external gravity (larger sq) requires a larger spatial 
scale for perturbations to become unstable. Note that 2/H^ = 2ttGT 1 q/c 2 s is equal to the Jeans 
wavenumber of a razor-thin disk (cf , Paper I) . 

We calculate Q c = k{2ttGp oc )~ 1 ^ 2 corresponding to marginal stability modes (uj = 0) from 




(18) 
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equation (18) and §3.1, and plot the results in Figure 36. Our simplified analytic approach gives 
Q c = 0.647 for unmagnetized disks with no external gravity (sq — ► oo), which is close to Q c = 0.676 
of Goldreich & Lynden-Bell (1965a) 's rigorous analysis. Critical values Q c for other sq and (3 are 
shown in Figure 36. We also perform axisymmetric (x—z) simulations using our numerical code to 
find Q c for realistic density distributions obtained in §3.1. The simulation results are displayed in 
Figure 36 as filled circles ((3 = oo), triangles (j3 = 10), and squares ((3 = 1), which are in excellent 
agreement with the results of the analytic approach. For average disk conditions (so = 1)> Qc 
is found to be ~0.75, 0.72, and 0.57 for (3 = oo, 10, and 1 cases, respectively. These represent 
approximately a 25% reduction relative to the corresponding Q c values of razor-thin disks. This 
suggests that the average disk with Q « 1.5 is quite stable to axisymmetric gravitational instability. 

We do not attempt in this paper to pinpoint threshold Qth values for gravitational runaways 
due to (nonaxisymmetric) swing amplification in finite-thickness disks, as we did in Paper I for 
razor-thin disks; there we found Q t h =1.3, 1.4, and 1.2 for (3 = oo, 10, and 1 cases, respectively. If 
we naively assume that Q c is lower by ~ 25% in nonaxisymmetric cases as well, then thick disks 
would have Qth ~ 0.9 — 1; our results from §4 would be consistent with such thresholds. The 
reduction of self-gravity in finite-thickness disks also implies that stable quasi-axisymmetric spiral 
shock solutions can exist at a stronger spiral arm potential than predicted in the razor-thin disks 
(cf, Paper II). 



4. Three-Dimensional Simulations 

To investigate local three-dimensional evolution of galactic gas disks under various conditions, 
we have performed a total of 10 numerical simulations. Table 1 lists the important parameters for 
each model. Column (1) labels each run. Column (2) gives the Toomre parameter Q (see eq. [8]). 
Note that Q < 1 in a thick disk does not necessarily mean axisymmetric gravitational instability 
(see §3 and Fig. 3). The relative strength of self-gravity compared to external gravity is given in 
terms of sq (see eq. [10]) in column (3); sq = 25 models strong spiral arm regions, while sq = 1 
characterizes average disk conditions. Column (4) lists the plasma parameter (3 characterizing the 
magnetic field strength (see eq. [9]). The shear parameter q (see eq. [1]) is given in column (5). 
Column (6) indicates whether self-gravity is implemented in the numerical evolution. Notice that 
even for the cases with "no" self-gravity for the evolution, the vertical total gravity profile from 
the "initial" equilibrium is included in order to keep the unperturbed equilibrium density profile 
intact. Columns (7) and (8) are our numerical resolution and the physical size of the simulation 
box in units of H. The relationship of H, Q, and so to physical values is given in §2.2. Finally, 
column (9) describes the relevant physics in each model. 
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4.1. The Parker Instability Without Self-Gravity 

In this section, we describe the nonlinear evolution of the Parker instability in the absence 
of self-gravity. These non-self-gravitating models demonstrate the effects of rotation and shear on 
the Parker instability, and will later be compared with the results of self-gravitating models. From 
the outset, we should mention that magnetorotational instability (MRI) in shearing disks with the 
present parameters is not important compared to the Parker instability. It is well known that MRI 
of toroidal magnetic fields grows much slower than MRI of poloidal fields (cf, Balbus & Hawley 
1998; Kim & Ostriker 2000). Even in a stratified disk with q = 1.5 and (3 = 12.5, toroidal MRI 
takes about ten orbital times to fully grow (Stone et al. 1996). Since the present models have 
weaker shear and stronger magnetic fields and are evolved only up to ~ (4 — 5)t or b, effects of MRI 
are negligible. 

4-1.1. The Parker Instability in Uniformly Rotating Disks 

We begin by considering the Parker instability under rigid-body rotation (model A). Figure 
4a shows the evolutionary history of the density maximum as a solid line. A three-dimensional 
visualization of the perturbed density with representative magnetic field lines is shown in the left 
panel of Figure 5. Snapshots of total surface density as well as x-y density slices are shown in 
Figure 6. 

Initially, density perturbations which are a superposition of numerous modes with different 
k adjust the disk by launching MHD waves. During this relaxation process, fluid motions arc 
very active at high-|z| regions where the density is low. Once the system finds the most unstable 
Parker mode, over-dense regions slide toward the midplane along the undulating field lines, while 
under-dense regions move upward at the expense of gravitational potential energy (cf, Hughes & 
Cattaneo 1987), causing the field lines to bend further. Using the evolution of vertical velocity 
dispersion, we find the growth rate of the most unstable mode is ~ 0Ac s /H in model A. After the 
maximum density p max = 1.75poo is reached at t/t rb ~ 1-8, the Parker instability stops growing 
and saturates; the corresponding maximum surface density is S max = 1.35£o- The primary reason 
for this is of course that the Parker instability is not a runaway process, but is eventually stabilized 
by tension forces from bent magnetic field lines (cf, Mouschovias 1974). 

As the left frame of Figure 5 and Figure 6a show, the dominant growing mode turns out to be 



4 Precisely speaking, the Parker instability that has the wavenumber vector k parallel to the initial magnetic field 
Bo is an undular mode of a more general class of magnetic buoyancy instabilities. The other branch of magnetic 
buoyancy instability is an interchange mode whose wavenumber vector is perpendicular to B (e.g., Hughes & 
Cattaneo 1987). Since the interchange mode requires 7 < 1 — 1/(2/3) to be unstable, where the z-independent plasma 
P is defined by eq. (9) and 7 is an adiabatic index (cf, J. Kim et al. 1998), our isothermal MHD models are all stable 
to the interchange mode. 
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antisymmetric across the midplane, and trailing with k x = 9(2ir/L x ) = 3.3/ H and k y = 2(2ir/L y ) = 
0.74/H; the corresponding wavelengths are X x = 2H and X y = 8.5H. The antisymmetric (in 
density) property of the dominant Parker mode in our model disks is consistent with the results of 
Santillan et al. (2000), who showed that numerical simulations with random initial perturbations 
generally exhibit antisymmetric modes, even though the linear growth rates for symmetric and 
antisymmetric modes of the Parker instability are almost indistinguishable in galactic environments 
(J. Kim et al. 2000). 5 

The dominant mode in model A is trailing (k x /k y > 0), as well. Linear dispersion relations 
for the Parker instability with rotation are independent of the sign of k x /k y (cf, Shu 1974), and we 
have found that for low-resolution simulations seeded with either leading or trailing modes, growth 
rates are similar. However, for the random perturbation form we adopt, trailing modes grow sooner 
because of a shorter readjustment phase, for the rotating cases. Additional tests for the nonrotating 
case show that the growth of trailing and leading modes is indistinguishable. 

Three-dimensional Parker instability attains its maximum growth rate at k x — > oo, because 
larger k x provides more space for a raised section of a flux tube to expand in the radial direction 
(Parker 1967). As pointed out by Shu (1974), however, the growth rate of k x = oo mode is not very 
different from that of the k x = modes and the growth rate is essentially constant for k x H £ 5. 
Since our high resolution simulations can resolve modes up to k x H ~ 12, the dominance of the 
k x H = 3.3 modes in model A is physical rather than numerical. The dominance of this particular 
mode over potential rivals may in part be attributed to the shape of initial perturbation spectrum 
which has larger power at smaller k x , while larger-/^ modes, if resolved, grow faster; a compromise 
may then select intermediate-/^ modes that have large enough initial power and fast enough growth 
rates. It is also possible that nonlinear effects help to select this mode, perhaps preferring close 
to the full disk thickness. We note that although uniform rotation tends to stabilize the Parker 
instability since the Coriolis force couples azimuthal and radial motions, its effect is not significant 
if kx/ky > 1 (cf, Shu 1974; Zweibel & Kulsrud 1975; Foglizzo & Tagger 1994). 

In addition to helping select trailing modes for initial growth, the effects of rotation are evident 
in other details of the model evolution. As material falls toward the midplane along wavy magnetic 
field lines, it also gains speed in the y direction, and its path is consequently deflected in the x 
direction by the Coriolis force. Gas moving toward a magnetic valley with positive (negative) v y is 
redirected radially outward (inward). The net result is that over-dense regions at valleys tend to 
rotate in the counterclockwise direction, which in turn twists magnetic field lines in the valleys (cf, 
Kim, Ryu, & Jones 2001). At the point when the instability begins to saturate, gas at the valleys 
is slightly overcompressed and begins to reexpand vertically. At the same time, the fields lines that 



5 Horiuchi et al. (1988) showed that the growth rates of antisymmetric and symmetric modes depend 
on their e parameter (see eq. [10] of Horiuchi et al. 1988). In galactic conditions, however, e « 
10 3 (7io/10kpc) 2 (r2/26kms _1 kpc _1 ) 2 (c s /7kms -1 ) 2 , so that their Fig. 5 reads that the maximum growth rates of 
antisymmetric modes are almost the same as those of symmetric modes. 
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previously became twisted due to the Coriolis forces now exert back torques, inducing gas motions 
in the clockwise direction. The meridional flows at the disk midplane, as displayed in Figure 6c, 
show this sort of behavior. The flows are generally subsonic near the midplane, while high altitude 
regions (\z\/H > 2.5) have supersonic velocities. 

Evolution of the system following saturation is rather complicated, involving nonlinear inter- 
actions of numerous MHD waves. Overpressured, dense regions shown in Figure 66 soon expand 
in the horizontal direction. Expanding material interacts with gas that moves along the vertical 
direction. These active, but still mostly subsonic, flow motions eventually homogenize the density, 
resulting in a smooth surface density distribution (Fig. 6d). The maximum surface density after 3 
orbits is only 1.13 times the initial value. 

The small-scale chaotic flow associated with nonlinear development of large- |k| modes is also 
likely to cause reconnection of the curved and twisted magnetic field lines, preferentially at low-|z| 
regions (cf, J. Kim et al. 1998; Kim, Ryu, & Jones 2001). This in turn reduces the field strength 
near the midplane, which was already decreased during the process of (midplane-crossing) Parker 
instability. 

We define the differential mass 



and plot dM/dz, d^s/dz, and dM/d^B as functions of z in Figure 7 for five time epochs. 
For t/^orb ~ 1) both profiles remain unchanged. As the Parker instability develops considerably 
(1 r$ t/t rb ~ 2), field lines float upwards, supplying magnetic flux at high altitude regions; mag- 
netic fields may escape the simulation box because of the outflow z-boundary conditions we adopted. 
Magnetic fields rise up and gas slips down, enhancing the mass-to-flux ratio further at the mid- 
plane (Fig. 7c), and evening out the z-distribution of magnetic flux (Fig. 7b), but mass enhance- 
ment near the midplane due to the Parker instability is generally quite moderate (Fig. 7a). At 
the end of the simulation, the overall density profile in model A is very close to another hydro- 
static equilibrium where the initial effective external force (g e ffi see eq. [17]) counterbalances the 
thermal pressure gradient alone. As Figure 7b suggests, the mean magnetic field (which is pre- 
dominantly azimuthal) is nearly constant with height, and therefore does not contribute much 
in the final force balance. At t/t rb = 4, the density- weighted velocity dispersions are calcu- 
lated to be a x = (pvD 1 / 2 /^) 1 / 2 « 2.5 x 10~ 2 c s , a y = (pv 2 ) 1 / 2 / (p) 1 / 2 « 7.3 x 10~ 2 c s , and 
o z = (/to 2 ) 1 / 2 / (p) 1 / 2 ~ 2.4 x 10~ 2 c s , where ( ) denotes a volume average. Given the low level 
of the velocity dispersion, the Parker instability is unlikely a significant source of interstellar tur- 
bulence in galactic disks. 




(19) 



and the differential magnetic flux averaged along the y-direction 




(20) 
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4-1-2. Effects of Differential Rotation 

One of the crucial impacts differential rotation makes on the evolution of the Parker instability 
is that it changes the power spectrum of initial perturbations. In a shearing disk, the wavefront 
of any disturbance evolves in response to the kinematics of the background flow; for uniform shear 
with q = — din f2/dln i?, the radial wavenumber of any pattern linearly increases with time as 
k x = k x (0) + qClkyt. In our chosen initial density perturbation spectrum, the mode with k x (0) = 
and k y = 1(2tt/L x ) = 0.37/H has much larger power than higher-/^ modes at the same k y . The 
higher-/^ modes, however, have higher growth rates than the k x = modes. As a consequence 
of the initial power distribution, non-shearing model A thus has to spend about one orbital time 
before the (large-A^) Parker modes grow enough to produce significant density fluctuation (Figure 
4a). On the other hand, the background shear in model B can transform the initial k x (0) = 
mode, without changing its power, into a higher-/^ disturbance - for which the Parker instability 
is more efficient. Therefore, in Figure 4a the initial growth of perturbations is more evident in 
model B than in model A. Structural analysis shows that this initially growing mode in model B 
is a z- symmetric (in density) mode with k y H = 0.37 (X y = 17 H) and confined to \z\/H < 1. In 
contrast, the most unstable mode in model A is a z- antisymmetric mode with k y H = 0.74. The 
growth rate of the k y H = 0.37 mode in model B is found to be ~ 0.2c s /H, which is just half of the 
growth rate of model A. 

While the initial evolution of shearing model B is governed by the growth of z-symmetric modes 
near the midplane, z-antisymmetric perturbations that have lower initial power (since k z is larger) 
also grow in high altitude regions at double the rate. As they grow, over-dense regions fall down 
toward the midplane and begin to interact with the symmetric modes. The right panel of Figure 
5 displays perturbed density and magnetic field lines in model B at t = 10fi -1 . The condensations 
evident at the midplane are associated with the symmetric mode with k x = W(2ir/L x ) (X x = 2H). 
Below or above these are "tadpole-shaped" condensations associated with the antisymmetric mode 
having k x = 20(2-7r/ L x ) (X x = H). Magnetic field lines from model B appear straighter than 
from model A, suggesting that the (longer-wavelength) symmetric modes partially prevent the 
(shorter-wavelength) antisymmetric modes from crossing the midplane. As antisymmetric modes 
grow dominant later on (t/i rb > 1-7), field lines become progressively more curved, but never to 
the extent seen in model A. 

For shearing model B, antisymmetric modes culminate their growth at t/i rb ~ 2.1. After 
model saturation, the system experiences vigorous interactions between symmetric and antisym- 
metric modes. An added complication is the background shear, which stretches out any conden- 
sation that forms. We plot snapshots of surface density X/So in Figure 8, which shows that at 
t/t or b = 2 the system is dominated by high-A^ modes, while surface density at i/t rb = 3 is rather 
uniformly distributed. To some extent, the smooth surface density is attributable to our limited nu- 
merical resolution. However, even if the true spatial scale at which dissipative processes smooth out 
strongly phase-mixed structures is small compared to our numerical resolution, there are evidently 
no nonlinear feedback processes to drive late-time growth of large-scale perturbations. Overall, the 
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late time character of the system is a new, large-scale vertical equilibrium containing small-scale, 
fluctuating density and velocity fields. The density-weighted dispersions of the perturbed velocity 
in model B at time t/t rb = 4 are a x ~ 0.10c s , <r y ~ 0.15c s , and a z ~ 0.07c s , which are too small 
to be of much significance for the observable interstellar turbulence in galactic disks. 

As mentioned before, symmetric modes in the shearing model B impede the midplane penetra- 
tion of antisymmetric modes, resulting in less inflated field lines. The antisymmetric modes have 
half the azimuthal wavelength of the symmetric modes, and also shorter radial wavelengths, which 
helps perturbations squeeze across the midplane. Nevertheless, the escape of magnetic flux into 
high altitude regions is significantly reduced. Figure 9 shows the variations of dM/dz, d^s/dz, 
and dM/d<&B for model B. Compared with Figure 7 for non-sheared model A, Figure 9 for model 
B illustrates that shear in the Parker instability results in a configuration where magnetic flux is 
enhanced in two layers (0.5H ^ \z\ ^ 1.5H) surrounding the midplane, while the corona has lower 
B (higher 0) than the non-sheared case (cf, Shibata, Tajima, & Matsumoto 1990). 6 Apart from the 
shearing background velocity, small-scale motions are chaotic, which can facilitate the reconnection 
of twisted magnetic field lines. The "layered" magnetic profile may in part arise from enhanced 
reconnection near the midplane, raising dM/d<&B at \z\/H ^ 0.5 (Fig. 9c). 

To summarize, we do not find any firm evidence that shear significantly stabilizes the Parker 
instability. This may be because Parker modes are the most unstable at large k x . By kinematically 
shifting the x-wavenumbers of perturbations at lowest k z and k y which have the largest initial 
amplitudes - and are associated with symmetric modes - the introduction of shear gives an early 
boost to symmetric modes. The more slowly-growing symmetric modes are later dominated by 
shorter-wavelength antisymmetric modes that grow more strongly (but start from lower ampli- 
tudes). The midplane-centered symmetric modes nevertheless interfere with the midplane-crossing 
antisymmetric modes, reducing the latter's amplification and limiting the inflation of field lines 
toward high latitudes. Except for having slightly more (less) magnetic support at moderate (low) 
\z\, the final state of the shearing model does not differ much overall from that of the nonshearing 
model, consisting of a new vertical equilibrium with enhanced central concentration. 

4.2. Gravitational Instability in Shearing Thick Disks 

4-2.1. Swing Amplification in Unmagnetized Disk 

Self-gravitating evolution of an unmagnetized, thick disk is generally very similar to that of 
a razor-thin disk, because without magnetic fields the induced motions are almost planar. One 
qualitative difference the disk's finite thickness makes is that it tends to be stabilizing by reducing 
self-gravity at the midplane (see §3). 



6 On the other hand, Miller & Stone (2000) showed that magnetorotational instability in a shearing disk tends to 
produce a strongly magnetized corona. 
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As shown in Figure 46 as a dotted line, the initial evolution of density in model C is governed 
by nonaxisymmetric swing amplification. With Q = 0.7 and sq = 25, hydrodynamic model C 
is stable to axisymmetric perturbations (Fig. 3). As wavefronts swing around from leading to 
trailing, epicyclic motion with the same rotational sense allows overdense fluid elements to linger in 
wave crests. This extended exposure to self-gravity enables perturbations to grow until the radial 
wavenumber becomes large (Goldreich & Lynden-Bell 1965b; Toomre 1981; see also Paper I). The 
swing amplification in model C saturates at i/i rb ~ 0.8, producing a state of nonlinearly interacting 
sheared wavelets of relatively high k x . Although the medium continues to shear kinematically, 
nonlinear interactions maintain a dominant population of wavelets at moderate k x . As displayed in 
the left panel of Figure 10, for example, the surface density map at t/t rb = 2 has the most power 
at k x = 7(2tt/L x ). 7 At this time, (|<fo y |) ~ 0.5(|v x |) and (\v z \) f» 0.04(|v x |), indicating the vertical 
motions are indeed insignificant compared with the planar motions. 

Much as we found for two-dimensional models (Paper I), the evolution of model C after satura- 
tion indicates that nonlinear interactions repopulate small-) k x \ modes, which then rejuvenate swing 
amplification. This second growth phase yields a few high-density filaments (middle panel of Fig. 
10). The filaments collide with each other and become gravitationally unstable, eventually forming 
two gravitationally bound clumps at the end of simulation. The mass of each clump is ~ 7% of the 
total mass, corresponding to ~ 4 x 10 6 M Q (right panel of Fig. 10); this is twice the corresponding 
two-dimensional Jeans mass, Mj = c^/(G 2 S ) = 2 x 10 6 M (c s /7.Okms~ 1 ) 4 (S o /65 Mqpc" 2 )- 1 , 
consistent with the results of Paper I where we found the fragments formed in a two-dimensional 
disk model with comparable Q and surface density typically have ~ (1/2 — 2)Mj (Model H06 in 
Paper I). 

In Paper I, we reported that when Q is sufficiently small, razor-thin disks form gravitationally 
bound clumps via three different types of secondary instabilities. These include parallel fragmen- 
tation (Q £j 0.8), gravitationally induced collisions of nonlinear sheared wavelets (0.8 £ Q £ 1-1), 
and rejuvenated swing amplification (Q £ 1.1 for unmagnetized disks). Had model C been evolved 
in the razor-thin approximation, it would have experienced gravitational runaway via parallel frag- 
mentation long before initial swing amplification saturated. It is remarkable that a thick disk with 
Q as small as 0.7 is barely unstable, requiring two successive phases of swing amplification before 
runaway occurs. 

4-2.2. Parker-Swing Instability 

The models of §4.1.2 and §4.2.1 represent systems where either gravity or magnetic fields are 
entirely turned off, to investigate the physics of sheared Parker and thick-disk swing instabilities 



7 Fourier analysis shows that modes of k x > W(2tt/L x ) have negligible power throughout the saturation stage, 
confirming that our simulation results are not limited by numerical resolution. 
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separately. Our primary interest here, however, is to understand how magnetic buoyancy may 
interact with the self-gravitating instabilities already known to exist from two-dimensional studies 
(Paper I). Model D focuses on the swing-Parker interaction, incorporating self-gravity, shear, and 
thermal-equipartition toroidal magnetic fields. Time evolution of the maximum density in model D 
is plotted in Figure 46 as a solid line. Clearly, the initial phase of evolution is dominated by swing 
amplification, with very similar behavior to model C (dotted line). Unlike in models A and B of 
the Parker instability without self-gravity, density in model D begins to grow immediately as k x 
increases. Embedded magnetic fields reduce the amplification factor in model D by 63% compared 
to model C. 

Swing amplification saturates at about t/t OI ^ ~ 0.8 when k x becomes sufficiently large. With 
embedded azimuthal magnetic fields, however, model D is also subject to the Parker instability. 
The shearing wavelets grown from the initial swing amplification in fact strengthen the symmetric 
Parker modes. Density in model D increases steadily for 1 i/i rb ~ 3; model B shows similar 
behavior at earlier times in Figure 4a. At the same time nonlinear interaction of the wavelets 
- as discussed in §4.2.1 - keeps the power spectrum centered at moderate k x . Figure 11 shows 
surface density snapshots at i/i or b = 2, 3, and 4.1 of model D. At i/i rb = 2, most of the power is 
concentrated in the mode with k x = 9(2ir/L x ) and k y = 1(2it/L x ). Vertical flows due to the Parker 
modes couple with in-plane (primarily radial) motions, enabling more vigorous nonlinear feedback. 
Rejuvenated swing amplification of small- \k x \ modes follows, and a few dense filaments form (middle 
of Fig. 11). They later experience mutual collisions, ultimately forming two gravitationally bound 
condensations, (right of Fig. 11). Each clump has 8% (left) and 12% (right) of the total mass, 
corresponding approximately to 9.6 x 10 6 M and 1.4 x 10 7 M Q . 

Because of the strong z-symmetric perturbations resulting from the initial swing amplification, 
model D does not provide fertile ground for growth of antisymmetric Parker modes. The latter 
may grow slightly in high- and middle-altitude regions, but never become dynamically important. 
Consequently, magnetic field lines do not float upward as much as in the non-self-gravitating model 
B. We draw dM/dz, d&s/dz, and dM/d<&B for model D in Figure 12, which illustrates that 
differential mass and flux at \z\/H ^, 1 are almost unchanged for i/t rb ~ 2. Nonlinear interaction 
of wavelets may cause magnetic field lines to reconnect, reducing d^s/dz near the midplane. 

Model D investigated how Parker instability interacts with swing amplification under condi- 
tions of relatively high gaseous surface density such that the (height-integrated) Q parameter is 
0.7 and gas gravity dominates stellar gravity locally (sq = 25). Regions in which both gas surface 
density and shear are so high may, however, be exceptional in contemporary galaxies. To inves- 
tigate the effects of the Parker-swing interaction under other galactic conditions, we also perform 
low-resolution simulations for two cases with lower gaseous surface density. Both of these models 
have so = 1; model G has Q = 1.5 and model I has Q = 1.0. For each of these, we also perform sim- 
ulations in which self-gravity is turned off (models F and H, respectively) to provide "pure Parker" 
comparison models. In addition, we turn off magnetic fields in model J to provide a Q = 1.0 "pure 
swing" counterpart to model I. The resulting evolutionary histories of maximum density in models 
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F-J are presented in Figure 13. 

When Q = 1.5 and sq = 1, corresponding approximately to average disk conditions, self-gravity 
is so weak that the growth of density via swing amplification in model G (solid line in Fig. 13a) is by 
less than a factor of two. Since the amplitude of the perturbed density field remains small, nonlinear 
feedback from wave interactions do not inhibit the kinematic increase of k x . The medium keeps 
shearing out, and Parker modes slowly emerge during the later part of evolution (i/i rb ~ 2.5). 
Model F, in which self-gravity is turned off (dotted line in Figure 13a), has late-stage evolution 
after saturation of Parker modes very similar to that of model G. Namely, a new quasi-static 
vertical equilibrium, with low-amplitude velocity and density fluctuations, is achieved. This state 
is similar to that of late-stage model B, except that the fluctuation level is lower by a factor of 
three. The reason for the differences between the two non-self-gravitating models (B and F) is that 
the former has a smaller scale height due to a stronger gas self-gravity in the initial equilibrium 
state 8 {HVL/c s curves in Fig. 1), thereby enabling the Parker instability to grow faster and more 
strongly (cf. eq. [12]). This is also why the clear transition from symmetric to antisymmetric Parker 
modes that occurs at t/t rb ~ 1-7 in model B (Fig. 4a, dotted line) is not obvious in the density 
history of model F (Fig. 13a, dotted line). The conclusive result from the Q = 1.5, sq = 1 model 
is that neither the Parker instability nor self-gravity could lead to significant density condensation 
under "average disk" conditions that are dominated by equipartition-amplitude, large-scale toroidal 
magnetic fields. 

Models H, I, and J investigate the Parker/swing connection under conditions of higher disk 
surface density such that Q = 1.0; we set sq = 1. The evolutionary histories of p mSLX /poo are shown 
in Figure 136. Evidently, comparison of the development of magnetized (solid line; model I) and 
unmagnetized (dashed line; model J) models shows that Parker modes can help push a region into 
self-gravitational runaway when Q is at a "marginal" level. Model I produces a couple of dense 
filaments at i/t rb = 4.8 and is expected to end in collapse, whereas the initial swing amplification in 
model J phase-mixes away within ~ 3 orbits. It is evidently a combination of the extra self-gravity 
produced in "valleys" of Parker modes combined with the vorticity-removing ability of magnetic 
fields that tips the balance in allowing a rejuvenated, runaway swing phase in model I. Model H, 
in which self-gravity is turned off (dotted line), shows the development of density perturbations for 
the "pure Parker" case, which is quite similar to that of model B (see Fig. 4a), even though sq and 
Q differ. 

In conclusion, the models of this section show that the Parker instability does not significantly 
alter the ability of a large-scale region of a galactic disk to form gravitationally bound clouds. Under 



8 As noted in the beginning of §4, our non-self-gravitating (for evolution) models still have the initial equilibrium 
self-gravity that together with external gravity balances thermal and magnetic pressure forces in the background 
state. Although the treatment of self-gravity is not self-consistent in model B, its primary function is as a control 
for the otherwise identical model D, in order to assess the small-scale evolutionary effects produced by self-gravity 
compared to Parker modes alone. 
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"average disk" conditions (Q = 1.5, sq = 1, q = 1), the disk turns out to be extremely stable, with 
only mildly fluctuating density in the final state. Both magnetized and unmagnetized disks can 
form bound clouds for the high gas surface-density conditions (Q = 0.7, so = 25) that represent 
spiral arm regions (except using "normal" shear). For intermediate surface density (Q = 1.0), 
only magnetized disks form bound clouds. While we have not pinpointed a threshold Q value 
for combined Parker-swing nonlinear runaway, models G (see Fig. 13a) and I (see Fig. 136) taken 
together bracket its value between 1.5 and 1.0 (when sq = 1). Even when bound clouds are created, 
the Parker instability plays only an auxiliary role by helping secondary swing amplification; the 
formation time scale and spacing of clouds are quite different from what the linear-theory Parker 
analysis predicts. 



4.3. Magneto-Jeans Instability 

The final issue we address in this paper is the development of the magneto-Jeans instability 
(MJI) in a three-dimensional disk. We previously studied MJI under uniform weak-shear conditions 
(Paper I), and under self-consistent spiral arm conditions (Paper II), both in a height-integrated 
two-dimensional idealization. While the latter varying shear and varying density approach is more 
realistic, in this first three-dimensional study we adopt the simplified uniform model, which is 
adequate to capture the basic physical processes: model E has the same high-gas-density-parameter 
set (Q = 0.7, so = 25) used in model D, but sets q = to represent better the reduced (or negative) 
shear that prevails within spiral arms (see Paper II). 

Physically, MJI is a form of nonaxisymmetric Jeans instability in rotating disks containing 
embedded azimuthal magnetic fields. Transverse (Alfvenic) magnetic perturbations produce radial 
tension forces that resist the stabilizing effect of the Coriolis force; magnetic transfer of angular 
momentum allows secular contraction instead of epicyclic motion (Lynden-Bell 1966; Elmegreen 
1987; Paper I). Although MJI is similar to the Parker instability in that both rely on background 
azimuthal magnetic fields and nonaxisymmetric perturbations, they are otherwise completely dif- 
ferent. The condition of low shear is essential for MJI, because magnetosonic waves begin to quench 
the Jeans instability when kinematic shear increases k x ; the Parker instability, in contrast, prefers 
high k x . Vertical motions are key in the Parker instability, whereas MJI involves primarily planar 
motions. Since both have been proposed as possible ways to trigger formation of giant clouds in 
spiral arms, it is interesting to investigate nonlinear evolution when both instabilities can occur. 

Based on model E, we conclude that MJI, if present, always dominates the Parker instability. 
Although equation (12) suggests the Parker instability can have a high initial growth rate, its 
development never produces highly nonlinear near-midplane density perturbations (see solid line 
in Fig. 4a). As the dashed line in Figure 46 shows, the density in model E (Q = 0.7) starts to grow 
immediately and this growth never saturates, just as we found for MJI in Paper I. We have also 
performed low-resolution simulations for MJI in nonshearing disks with Q = 1.0 and 1.5 (not listed 
in Table 1), and confirmed that MJI can drive them into gravitational runway. The time scale for 
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bound cloud formation becomes longer as Q increases, so that disks with Q = 1.0 and 1.5 form 
dense clumps within 1 and 2 orbits, respectively, while model E with Q = 0.7 takes less than one 
orbital time to become fully unstable (Fig. 46). 

In Figure 14, we plot perspective views of isodensity surfaces and magnetic field lines, together 
with midplane distributions of density, velocity, and magnetic fields, for model E at t/t rb = 0.5. 
As the perspective view shows, the magnetic fields run almost parallel to the x-y plane near z = 0, 
as do the velocity vectors even during highly nonlinear stages. Centered at (x/L x ,y/L y , z/L z ) = 
(—0.10, —0.19,0) and (0.14,0.05,0), two condensations are produced within one orbital time; these 
each contain 27% and 20% of the total mass, corresponding roughly to ~ 3 x 10 7 M Q . 

As discussed in §3.2 for axisymmetric gravitational instabilities, the dilution of in-plane self- 
gravity due to the finite scale height also affects the most unstable scale lengths and growth rates 
of MJI. The simple replacement Sofc — ► 2pooHk/(l + Hk) is sufficient to take into account this 
geometrical dilution effect in an approximate fashion. Utilizing equation (23) of Paper I (see also 
Lynden-Bell 1966), the instantaneous dispersion relation for MJI in thick disks becomes 
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k x (0) + qQk y t. Perturbations require k y ^ (nonaxisymmetric) 

1)/2H to be unstable to MJI (see eq. 



where k 2 = k x (t) 2 + k 2 and k x (t) 

and £#00(1 + kH) < 2, or k < k c = (a/1 + 8H/H 
[15] for the definition of Hoc). Note that the critical wavenumber k c of MJI is independent of 
magnetic pressure (except through H) and rotation, so that MJI operates at smaller scales along 
the magnetic field than axisymmetric gravitational instability. Since H/H^ ~ 1 for so ~ 1 (Fig. 
1), the condition for MJI is roughly kH ^ 1, which is similar to the instability criterion for the 
Parker modes. 9 

Using equation (21) and the results of §3.1, we obtain the peak (k x = 0) growth rates for MJI 
for model E (Q = 0.7, sq = 25 and (3 = 1) as a function of k y , and plot them in Figure 15. For 
comparison, we also plot the growth rates for the corresponding razor-thin disk. The maximum 
growth rate for a thick disk is ~ 1.12f2 at k y>iaax = 0.74/ H, so Xy !iaax ~ 8H. Since L y ~ 17 H for our 
box, this would predict \ y , ma , x /L y ~ 0.5, i.e., two x-aligned filaments. This is indeed consistent with 
the density structure in Fig. 14, although there is also a small contribution from the A^max = L y 
mode that has higher power from our initial conditions. Note that the most unstable A^ of MJI 
is in fact the same as that of the dominant Parker mode that grows in non-sheared model A, 
although that mode is midplane-antisymmetric in density. For both x- and y-scales comparable to 



Aj/,max, the corresponding mass is M = So(27r /k y>1 



3 x lO 7 M (E o /65 M pc- 2 )(i7/8Opc) 



9 In investigating the formation of filamentary structures in the ISM, Chou et al. (2000) studied disk models that 
are similar to model E, that is, magnetized, uniformly-rotating, and self-gravitating. They identified the structure 
forming mechanism with the Parker- Jeans instability, but we believe that while the Parker mechanism may contribute 
during the linear growth phase, the primary driver of nonlinear structure formation is the MJI. 
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consistent with the masses of clumps formed in model E. From Figure 15, we note that a razor-thin 
disk would characteristically respond to MJI twice as fast, at half the length scale, of a thick disk. 
Consequently, the characteristic Jeans mass is about four times bigger in a thick disk than the thin- 
disk prediction. We note that this result is consistent with simple prediction for an unmagnetized 
disk using the vertical equilibrium condition of equation (15): the central density is 7tGSq/2c^, and 
with the three-dimensional Jeans wavelength Aj^d = c s (ir jGp) 1 ! 2 , the nominal three-dimensional 
Jeans mass would be Mj^d = 7r\/2cf /(G 2 S ). With Aj^d = c^/GHq, the two-dimensional Jeans 
mass is Mj i2D = 4/(G 2 S ) ~ 0.23M Ji3D . 

5. Conclusions 
5.1. Summary 

A variety of evidence supports the idea that giant, star-forming clouds form as the result of 
large-scale collective effects - i.e., dynamical instabilities. Candidate mechanisms that have been 
proposed on the basis of linear-theory analysis include the swing amplifier (Goldreich & Lynden- 
Bell 1965b; Julian & Toomre 1966), magneto- Jeans instabilities (Lynden-Bell 1966; Elmegreen 
1987; Paper I), and the Parker instability (Parker 1966). The first two of these involve primarily 
in-plane motions and are fundamentally driven by self-gravity, with either shear or magnetic tension 
acting to neutralize the stabilizing tendency of the Coriolis force. For the third mechanism, vertical 
motion, driven by magnetic buoyancy, is crucial. 

Based on linear-theory growth rates and characteristic spatial scales, all of the above mecha- 
nisms could in principle be dynamically important under various galactic conditions. In practice, 
because giant clouds represent highly overdense regions compared to mean ISM properties, it is im- 
portant to understand how the candidate instability mechanisms develop in the nonlinear regime: 
i.e., whether and/or how an instability saturates. Addressing these questions requires direct nu- 
merical simulations to evolve the fluid dynamics equations into the nonlinear regime. 

In previous work, we have investigated the nonlinear development of the swing and magneto- 
Jeans mechanisms in a thin-disk (two-dimensional) approximation (Papers I, II). Here, we extend 
those studies by allowing for fully three-dimensional dynamics, which also enables us to investigate 
the nonlinear development of the Parker instability allowing for rotation, shear, and self-gravity. 
Our primary goals were to compare dynamical development, and especially "final-state" outcomes, 
of model disks under various conditions. By controlling for different physical effects (i.e., turning 
gravity, rotation, shear, and magnetic fields "off" and "on"), and exploring parameter space, we 
were able to assess the potential consequences for bound cloud formation of the various instabilities 
alone and in combination. 



Because H varies inversely with E from eq. (15), the condensation mass would be proportional to E 
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The numerical models we adopt are radially localized, vertically stratified, shearing, isothermal 
disks with embedded magnetic field lines (see §2.1). The initial rate of shear in the background 
flow is measured by q = — <iln Q/d\n R, and this overall shear gradient is maintained by enforcing 
sheared-periodic radial boundary conditions; we study cases with q = 1 and 0. The initial magnetic 
field is assumed to be azimuthal and its vertical stratification is determined such that the Alfven 
speed va is spatially constant. The density and magnetic field strength in our model disks are 
characterized by the Toomre parameter Q (see eq. [8]) and the plasma parameter f3 (see eq. [9]), 
respectively. Magnetic fields are either initiated at an equipartition level (va/c s = l) 11 , or set to 
zero. We set up initial magnetohydrostatic equilibria with thermal and magnetic pressure gradients 
balancing gaseous self-gravity and external gravity from a stellar component. We adopt a linear 
model for the external (stellar-disk) gravity (see eq. [7]) and parameterize its magnitude using so 
(see eq. [10]), which corresponds roughly to the ratio of the gas disk's vertical gravity to the external 
gravity at one scale height. The fiducial sets of parameters our model disks adopt are Q = 1.5, 
so = 1 for the "average disk", and Q = 0.7, so = 25 for highly compressed spiral arm regions; we 
also perform additional simulations with Q = 1 and so = 1. 

Our simulation boxes are 8 scale-heights (H) thick, with area either 17H x 17H or 25H x 25H 
in the x—y plane. The relation between dimensionless simulation parameters and physical scales 
is explained in §2.2. For all our simulations, we introduce a spectrum of low-amplitude random 
perturbations, and numerically integrate the ideal MHD equations up to four or five orbital periods. 
In §3, we first constructed the density profiles for vertical static equilibria and calculated the 
variation of disk thickness as a function of So (§3.1), and then studied axisymmetric gravitational 
instability of these thick disk models (§3.2). Finite disk thickness reduces the critical value of the 
Toomre stability parameter relative to the zero-thickness limit, and also increases the spatial scale 
required for instability. We found that for so = 1, Q c ~ 0.75, 0.72, and 0.57 for (3 = oo, 10, and 1 
cases, respectively, suggesting that the average disk is highly stable to axisymmetric perturbations. 
We also showed that a simple modification of the axisymmetric dispersion relation to account for 
finite thickness (eq. [18]) yields excellent agreement with the results of simulations for the value of 

Qc 

In §4, we present nonlinear evolution of fully three-dimensional disk models, first considering 
the Parker instability in isolation (§4.1), then investigating evolution of the swing and Parker-swing 
mechanisms (§4.2), and finally considering how self-gravitating perturbations grow under low shear 
conditions (§4.3). 

The main results drawn from these three-dimensional numerical simulations can be summarized 
as follows: 

1. In the absence of self-gravity, the dominant growing Parker mode under uniform rotation 
is found to be antisymmetric in density with respect to z = and tailing with radial wavelength 



Thus magnetorotational instability cannot grow; see below. 
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X x = 1.9H and azimuthal wavelength X y = 8.5H, and has a growth rate ~ 0Ac s /H. Density fluc- 
tuations grow as over-dense regions slide toward the midplane along undulating magnetic field lines 
and unloaded portions of the field buoyantly rise, but growth slows as nonlinear saturation sets in. 
The maximum density at saturation is less than twice the initial midplane density. After saturation 
of the dominant Parker mode, nonlinearly developed small-scale modes interact with each other, 
dispersing density structures characteristic of the Parker instability. Escape of magnetic flux across 
the vertical boundaries and likely reconnection of the twisted magnetic fields prompted by small- 
scale motions raise the mass-to-flux ratio near the midplane. The system progressively moves to a 
new global vertical equilibrium in which the magnetic field is nearly uniform, hence magnetic pres- 
sure support against gravity is negligible. In the saturated state, small scale fluctuating velocities 
are still present, but their density-weighted dispersions are less than 10% of the isothermal speed 
speed. 

2. When a disk is subject to differential rotation, midplane symmetric (in density) Parker 
modes with X y = 17H are the first to dominate the evolution. This initial dominance is explained 
by the large initial amplitude of longer-wavelength perturbations, combined with the ability of 
shear to increase k x to the regime where Parker instabilities are more efficient. Smaller-Ay anti- 
symmetric modes start with lower amplitudes than the symmetric modes but grow more rapidly, 
moving gradually downwards from high altitudes to interact with the symmetric mode. The mode 
interaction reduces magnetic field line inflation and loss compared to the non-shearing case. The 
final state of the non-self-gravitating shearing disk model is similar to that of the non-shearing disk, 
consisting of an overall vertical equilibrium with small-scale, low-amplitude fluctuations in density 
and velocity fields. 

3. When we include both self-gravity and uniform shear in our models, the nonlinear evolution 
of structure is principally dominated by the swing amplification mechanism, although Parker modes 
play a role in seeding structure. Low-Q (Q = 0.7) disk models end in gravitational runaway, while 
high-Q (Q = 1.5) disk models remain extremely stable. In both high and low Q limits, the ultimate 
outcome is decided independent of the presence or absence of magnetic fields. When Q is marginal 
(Q = 1), on the other hand, we find that only magnetized models form bound condensations, 
suggesting that Parker instability may play a supplementary role in destabilization. The Parker 
instability raises the level of density fluctuation and may help channel power in the saturated 
state into small wavenumber modes. It is the "rejuvenated" swing amplification (see Paper I) that 
ultimately drives the system to form dense condensations, however. For cases in which bound 
clouds do form, the sizes and spacings are not characteristic of Parker modes, but comparable to 
the results from thin disk models (with typical masses M ~ (2 — 3) times the thin disk Jeans mass). 

4. Weakly shearing thick disks, like thin disks, are susceptible to the magneto-Jeans instability 
(MJI), in which tension forces from azimuthal magnetic fields vitiate the stabilizing Coriolis force. 
The azimuthal wavelength criterion for MJI in finite-thickness disks is similar to that of the Parker 
instability (k y H 1), and modal growth rates are also similar. Our nonlinear simulations show, 
however, that MJI under weak-shear (spiral arm) conditions is ultimately much more violent than 
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the Parker instability, and can form dense condensations within one orbital time. This is because 
the Parker instability is self-limiting, whereas the MJI is a runaway process. The mass of dense 
condensations amounts typically to the thick-disk Jeans mass (~ 4 times the thin-disk Jeans mass). 

5.2. Discussion 

The numerical simulations in this paper suggest that formation of bound clouds cannot occur 
as a direct consequence of Parker instability, although Parker instability may play a role in seeding 
swing amplification in high-shear environments or MJI in regions of weak shear. There are several 
reasons why the Parker instability is ineffective at producing large, dense clouds. First, although 
the growth time of the Parker instability from linear-theory analyses is a few times 10 7 yrs (cf, Kim 
& Hong 1998), the process does not proceed in a runaway fashion, but is stabilized at a relatively 
modest density by magnetic tension forces (e.g., Mouschovias 1974). Second, while the undular 
Parker instability prefers wavelengths in the azimuthal direction X y ~ 2ttH (i.e., k y H ~ 1 in rough 
agreement with observed spacing between molecular complexes) , there is no similar preferred scale 
in the radial direction (e.g., Parker 1967; Giz & Shu 1993; J. Kim et al. 1998). Growth rates increase 
slowly with k x , and are nearly uniform for k x H ^ 5. Indeed, even though the initial perturbation 
amplitudes in our simulations increase toward large scale, we find that small- X x modes dominate 
the evolution. The results is long filamentary structures with a radial spacing about 10 times 
smaller than the azimuthal spacing in surface density maps during the initial growing phase (see 
Fig. 6). 12 Because the perturbed self-gravitational potential 5<fr s oc —Sp/k 2 is increasingly small 
for large k x , effects of self-gravity on development of the Parker- Jeans instability are weak for the 
dominant nonlinear modes. 

Third, while two-dimensional simulations for the Parker instability inarguably reach Mouschovias 
(1974)-type final static equilibrium configurations (cf, Santillan et al. 2000), three-dimensional sim- 
ulations (J. Kim et al. 1998, 2000, and this work) indicate that such equilibrium configurations are 
transient at best, evolving rapidly into a state where density is relatively smooth in the horizontal 
direction. It is unclear whether this transition is initiated by an instability of the static undulating 
configurations, as claimed by Lachieze-Rey et al. (1980), by active magnetic reconnection near the 
midplane caused by small-scale chaotic motions, or simply by comparable nonlinear saturation lev- 
els of many high-fc^ modes. When shear is present, kinematically driven phase mixing contributes 
to washing out large-scale radial structures. In any case, the three-dimensional simulations show 
that it is almost impossible to find any signature of the undular Parker instability in the later stages 
of evolution. 

The result that Parker instability cannot directly account for GMC formation has previously 



The corresponding radial wavelengths are one or two times the disk scale height for symmetric and antisymmetric 
structures, respectively, possibly indicating nonlinear selection of these modes. 
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been emphasized by other authors. Kim, Ryu, & Jones (2001), for example, find only a maximum 
factor of 1.5 enhancement in surface density in their three-dimensional non-self-gravitating simu- 
lations with solid-body rotation, and conclude that Parker mode alone cannot be responsible for 
forming molecular clouds (see also e.g., Basu, Mouschovias, & Paleologou 1997; J. Kim et al. 1998, 
2000; Santillan et al. 2000). These authors speculate on the possibility of the Parker instability 
coupled with other physical processes such as self-gravity, differential rotation, radiative cooling, 
and spiral arm potentials. We show in the present paper that self-gravity and galactic differential 
rotation are not much help. Therefore, we argue that the Parker instability, even combined with 
gaseous self-gravity and differential rotation, cannot be regarded as the primary formation mech- 
anism for giant molecular clouds at least in galactic disks at large. While effects on the Parker 
instability of stronger self-gravity and rapidly varying shear inside spiral arms have yet to be stud- 
ied in detail, we do not anticipate significant differences in the main conclusions. Although the 
width of a spiral arm might help in promoting Parker modes of similar X x over shorter-wavelength 
competitors (cf., Franco et al. 2002), it is not yet clear how effective this can be if the arm transit 
time is comparable to the growth time. 

While larger or smaller shear rates do not seem to alter the basic character of Parker instabilities 
appreciably, the shear rate is crucial for selecting which type of (in-plane) gravitating amplifier can 
operate. In a weakly shearing regions, swing is absent, but the MJI is very powerful at initiating self- 
gravitational condensations. In Paper I, we showed MJI is a potential means to prompt observed 
active star formation near the central parts of galaxies where rotation curves are nearly rigid- 
body (see also Elmegreen 1987). In Paper II, we argued that observed spur structures jutting 
radially outwards from spiral arms can be understood as a direct consequence of the MJI occurring 
within the arms. A spiral arm potential has three important effects on the gaseous medium. 
(1) By compressing gas, the Jeans length - which determines the preferred MJI scale along the 
mean magnetic field - is reduced; the growth time for self-gravitating modes also drops. (2) The 
spatially-varying density distribution causes a gradient in the shear profile perpendicular to the 
arm; since transiting gas first experiences reversed shear and then returns to normal shear, the 
spatially-varying profile tends to keep radial wavenumbers small within the arm, which enhances 
the efficiency of MJI. (3) It reduces the radial scale length for mass collection by MJI to at most 
the arm width, so that the masses of clumps produced are limited. 

Compared to our previous results from razor-thin disks, the weaker self-gravity in three di- 
mensions could imply somewhat more massive (~ 10 7 M ) clouds would tend to form from MJI 
under realistic conditions. On the other hand, the weaker self-gravity may allow stronger spiral 
potential perturbations that give a higher density contrast between arm and interarm regions in 
stable quasi-axisymmetric spiral shock profiles (cf. Paper II). The Jeans mass at the spiral arm 
density peak could be correspondingly smaller, producing lower mass clouds in three-dimensional 
spiral arm regions. Since the marginal wavenumber of MJI parallel to B is about the same as that 
of the Parker instability in spiral arm conditions (see §4.3), one might expect some cooperation 
between MJI and the Parker instability inside spiral arms. The present models under uniform 
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in-plane background conditions do not point to significant modifications of the nonlinear MJI by- 
Parker modes. Effects from Parker-MJI interactions could be more prominent with a realistic spiral 
arm model (cf. Martos & Cox 1998), however, and it would be particularly interesting to learn 
whether the final mass scale changes markedly compared to predictions from uniform disk models. 

By focusing on cases with relatively strong (/3=1) equilibrium azimuthal magnetic fields in 
our MHD models, we have concentrated solely on effects of the Parker instability and MJI in 
the present work. An additional dynamical process prevalent in magnetized disk systems is the 
magnetorotational instability (MRI; Balbus & Hawley 1991, 1998). Because MRIs are most vigorous 
in their small-scale growth when the thermal pressure is large compared to magnetic pressure 
(/3 > 1), they have been studied primarily in the context of accretion disks. However, Sellwood 
and Balbus (1999) have argued that the MRI may also be important in producing turbulence 
in galactic disks where the effective value of (3 is close to unity, citing as evidence the significant 
nonthermal velocity dispersions (~ 6 km s _1 ) present in outer H I disks with low local star formation 
rates and apparent gravitational stability. It is of great interest to investigate whether coupling 
between the MRI and self-gravitating modes under galactic conditions could represent an important 
new mechanism for forming large-scale ISM condensations. 

We are grateful to an anonymous referee and N. Turner for thoughtful comments and to S. 
White, T. Sano, and J. Lee for help in using IDL to make three-dimensional visualization. This 
work was supported by NASA grant NAG 59167. 
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Table 1. Parameters of Three-Dimensional Simulations. 
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Note. — See §4 for explanation. 
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Fig. 1. — Variation of the vertical scale height H of the gas distribution in a magnetohydrostatic 
equilibrium with varying gaseous-to-stellar surface density (oc s^ 2 ). The H/Hoo curves show the 
thickness relative to that of a purely self-gravitating gas disk, while the HQ/c s curves show the 
thickness relative to a fixed fiducial value. Strong self- or external gravity leads to a smaller scale 
height. Magnetic pressure makes a disk thicker. 

Fig. 2. — Height dependence of driving force g e R (solid lines) of the Parker instability for two 
model galactic disks. Dotted and short-dashed lines correspond to self-gravity and external gravity, 
respectively, (a) Under average disk conditions (Q = 1.5 and sq = 1), gaseous self-gravity is 
comparable to the external gravity up to z/H ~ 1, beyond which the latter dominates. For 
comparison, we plot the external gravity models adopted by Giz & Shu (1993) and Kim & Hong 
(1998) using dot-dashed and long-dashed lines, respectively, (b) For strongly-compressed gas in a 
spiral arm (Q = 0.7 and So = 25), self-gravity is much more important than external gravity. 

Fig. 3. — (a) The most unstable wavenumber fc max and (b) critical Q c for axisymmetric gravitational 
instability in a finite-thickness disk. Various symbols (circles for (3 = oo, triangles for (3 = 10, and 
squares for (3 = 1) plotted in (b) are obtained by numerical simulations, showing excellent agreement 
with the analytic estimates. See text for details. 

Fig. 4. — Time evolution of maximum density in models A - E with Q = 0.7 and sq = 25. (a) Models 
without self-gravity show only moderate increase of density, never developing highly nonlinear 
perturbations, (b) With self-gravity included, the disk models eventually become gravitationally 
unstable, forming bound clumps at the end. 

Fig. 5. — Perturbed density (color scale) in three different planes (x = —0.5L X , y = —0.5L y , and 
z = 0.12L Z ) and selected magnetic field lines for (left) model A and (right) model B at t/t rb = 1.6. 
Only the portions of magnetic fields that start from x = —0.5L X and z = 0.12L Z and rise above 
the z = 0.12L Z plane are shown. The midplane (z = 0) is marked by thin line in each box and 
colorbars label 5p in units of the initial central density. 

Fig. 6. — Snapshots from model A of the Parker instability with uniform rotation, (a) Total surface 
density map (x—y plane) at i/t rb = 2. (b) Density structure at the midplane (z = 0) at t/t rb = 2. 
(c) Square section in (b) is enlarged to show density and horizontal velocity vectors in detail. The 
arrow outside the box shows the sound speed for comparison, (d) Surface density map at t/t rb = 3. 
Gray scales are linear in units of So or poo- 

Fig. 7. — (a) Differential mass dM/dz and (b) magnetic flux d&s/dz, averaged aver the y-axis, 
and (c) differential mass-to-flux ratio dM/d&B of the Parker instability simulation with uniform 
rotation (model A) are drawn as functions of the vertical height z. 
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Fig. 8. — Surface density maps projected on the x-y plane of the simulation for the Parker instability 
including shear (model B) at i/i rb = 1) 2, and 3 from left to right. Gray scale bars label S/Sq- 



Fig. 9. — Same as Figure 7 except for model B that contains shear. 



Fig. 10. — Snapshots of surface density of the hydrodynamic model C at t/t or \, =2.0, 2.9, and 3.4 
from left to right. Numbers labeling gray-scale bars correspond to logS/So- 



Fig. 11. — Snapshots of surface density of model D at t/t OT ^ =2.0, 3.3, and 4.1 from left to right. 
Numbers labeling gray-scale bars correspond to logS/So- 



Fig. 12. — (a) Differential mass dM/dz and (b) magnetic flux d^s/dz, averaged aver the y-axis, 
and (c) differential mass-to-flux ratio dM/d^B of model D. 



Fig. 13. — Time evolution of maximum density for (a) models F and G with Q = 1.5 and (6) models 
H - J with Q = 1.0. In this plot, G^O indicates self-gravitating cases, while non-self- gravitating 
cases are marked by G = 0. 



Fig. 14. — Left: Perspective visualization of isodensity surfaces (p/poo = 1-5) and selected magnetic 
field lines in model E at t/i rb = 0.5. Starting positions of magnetic field lines lie at y/L y = —0.5 
and z/L z = —0.45, 0, and 0.45 from bottom to top, respectively. Right: Density in logarithmic 
gray scale, velocity vectors (arrows), and magnetic field lines (thick solid lines) at the midplane in 
model E at t/t rb = 0.5. Note that z-components of velocity and magnetic fields are almost zero 
near the midplane. 



Fig. 15. — Growth rates of magneto-Jeans instability for thick and thin disk models with Q = 0.7, 
so = 25, (3 = 1, and k x = 0. Overestimated self-gravity in thin disks implies larger growth rates 
and smaller unstable length scales. 
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